Comparison of Kinetic Models Applied for Transport Description in Polymer Inclusion Membranes

Five mathematical models for transport description in polymer inclusion membranes (PIMs) were presented and compared via regression analysis. The applicability of the models was estimated through the examination of experimental data of Zn(II), Cd(II), Pb(II), and Cu(II) ions transported by typical carriers. In four kinetic models, a change in the feed and stripping solution volume was taken into account. The goodness of fit was compared using the standard error of the regression, Akaike information criterion (AIC), Bayesian (Schwarz) information criterion (BIC), and Hannan–Quinn information criterion (HQC). The randomness distribution in the data was confirmed via a nonparametric runs test. Based on these quantities, appropriate models were selected.


Introduction
In past decades, polymer inclusion membranes (PIMs) have been frequently used in analytical applications. PIMs are employed to separate and preconcentrate various species such as metals, organics, and pharmaceutical compounds [1][2][3][4][5]. Their use in analytics requires the development of quantitative descriptions of the transport kinetics, which can be obtained through mathematical modeling. Appropriate equations can be used for the calculation of time-dependent values of, for example, the selectivity, concentration, or recovery factors, since an accurate prediction of these quantities is important in analytics, especially in the procedures of sample preparation [6].
For many years, simple mathematical models similar to those describing the reaction kinetics have been used to describe transport in many membrane systems. These equations have been primarily applied for describing transport in liquid membrane (LM) and PIM systems. As their main advantage, they enable the calculation of the parameters that influence transport efficiency (e.g., maximum flux) and predict concentration changes in respective solutions of the adjacent membrane solutions.

Model No. 1
Despite the availability of many different kinetic equations, one of the most commonly employed mathematical models is that proposed by Danesi [7], whose mathematical form is identical to the differential equations describing the kinetics of first-order reactions: and results from the general transport scheme in the form of In the above equations, A denotes the membrane surface area [cm 2 ]. This simple model, as represented by Equations (1) and (2), is most frequently used to describe transport kinetics in PIMs. For example, only last year, it was used to describe the removal of fluoride [8] and phenol [9], the separation of mercury(II) [10], the separation of lithium and magnesium [11,12], the recovery of bismuth(III) [13] and scandium [14], the extraction of arsenic(V) [15], the separation of Pb(II), Zn(II), and Cd(II) ions [16], and the removal of antibiotics [17].
Noteworthy is the fact that, in some cases, the model fit quality is poor [15,17] or does not satisfy the criterion of a random distribution of residuals [10,16]. This means that a different or more advanced model should be used to describe the transport kinetics.

Model No. 2
It was previously shown that the transport of Cd(II), Zn(II), Pb(II), and Cu(II) ions through PIMs containing reactive ionic liquids or D2EHPA (di-(2-ethylhexyl)phosphoric acid) as a carrier obeys the kinetic laws of a reversible first-order reaction [18,19]: According to the general transport scheme above, the rate of change in concentrations in the respective solutions during the membrane transport can be expressed by This model was also successfully applied for the description of Pb(II) ion transport through PIMs containing calixresorcin [4] arene derivatives as carriers [20].

Model No. 3
In the models presented above, the initial maximum flux can be interpreted as a flux related to the sorption of molecules into the membrane. In this case, the value of the initial flux into the receiving phase has no physical meaning because in membrane processes controlled by diffusion, a time lag should be expected. This means that the minimum period of time is needed for the formation of the complex and its penetration into the membrane interior until the release of the first molecules into the receiving phase. At the beginning of the process, the zero value of the transported substance flux into the receiving solution should be expected if the membrane has not been pre-equilibrated with the feed solution. The application of effective carriers and the appropriate membrane composition allow for minimizing the time lag effect to such an extent that it is not observable in the process. However, in order to describe transport through PIMs, considering the diffusive nature of the process, the time lag should be taken into account, as it is always present in this type of system. For this purpose, a model whose equations are similar to the equations describing the kinetics of first-order consecutive reactions can be used: Such a mathematical model is mainly applied for the description of the transport of various substances through bulk liquid membranes (BLMs). For instance, this model was applied in the case of the pertraction of chromium(VI) [21,22], chromium(III) [23], chromium(II) [24], mercury(II) [25], cobalt(II) [26], cadmium (II) [27], and cyanide and thiocyanate ions [28,29].

Model No. 4
Analogous to the model proposed by Danesi [7], in this model, a problem with its fitting to experimental data, as well as the related incorrect residual distribution, is observed [21,23,28]. The quality of the model fit can be improved by applying a model similar to a consecutive reaction scheme with a slow reversible step, antecedent to the irreversible step, according to the scheme which can be described by the following set of differential equations: This model was previously applied for the description of phenol [30], L-isoleucine [31], and strontium(II) pertraction [32] through BLMs.
In summary, there are at least four simple kinetic models that can be used to describe transport through PIMs. Nevertheless, only the model proposed by Danesi [7] is most frequently applied.

Model No. 5
It should be noted that in the case of the linear concentration vs. time dependence for transported substances, an equation similar to those describing a zero-order reaction (which corresponds to stationary or pseudo-stationary conditions of the transport) can be used: Despite the many different kinetic models that enable the determination of the initial maximum flux (J M ) and a quantitative description of the concentration change in the feed and receiving solutions, some authors calculate J M from the first derivative of an exponential decay function fitted to the feed solution concentration dependence [33][34][35][36][37]. This approach probably results from problems occurring when fitting the most commonly used kinetic model (model No. 2) to the experimental data. The calculated J M value describes the system efficiency. However, it is impossible to obtain a quantitative description of time-dependent concentrations, especially in the receiving solution. In this case, there is also no possibility of a physicochemical interpretation of the exponential decay function, contrary to the typical kinetic models. The main aim of this work is, therefore, to discuss the applicability of the above simple kinetic models and to select the most appropriate model for the transport description of substances through PIMs. The usefulness of the models was evaluated with the use of typical parameters describing the quality of the model fit to the experimental results of Zn(II), Cd(II), Pb(II), and Cu(II) ion transport. In the model calculations, a change in the feed and stripping solution volume (because of sampling) was taken into account. The goodness-of-fit evaluation was estimated using the standard error of the regression, Akaike information criterion (AIC), Bayesian (Schwarz) information criterion (BIC), and Hannan-Quinn information criterion (HQC). A nonparametric runs test was also used to examine the randomness of the residuals.
The novelty of this research is the development of a method for the selection of a proper model based on appropriate fit quality parameters and a runs test as well as the application of kinetic models that have never been used for the description of transport through PIMs.

Experimental
Experimental studies of Zn, Cd, Pb, and Cu ion transport through PIMs containing various types of carriers were carried out. The following carriers were used for this purpose: TOPO (tri-n-octyl phosphine oxide, 90% Sigma Aldrich, St. Louis, MO, USA), D2EHPA (di-(2-ethylhexyl) phosphoric acid, 97% Aldrich, St. Louis, USA), Aliquat 336 (methyl trioctyl ammonium chloride, Aldrich St. Louis, MO, USA), Cyphos IL 101 (trihexyl(tetradecyl)phosphonium chloride, >97% Solvionic, Toulouse, France), and RILC8_Br (3-(1,3-diethoxy-1,3-dioxopropan-2-yl)-1-octylimidazolium bromide, synthesis described in [18]). Experimental studies were carried out in the system described in detail in [18,19]. The feed solution was composed of Zn(II), Cd(II), Cu(II), and Pb(II) nitrates (Sigma-Aldrich, St. Louis, MO, USA, reagent grade, purity ≥ 98%) dissolved in 0.5 M HCl (200 cm 3 ) with an initial concentration equal to 0.002 M. Only in the system with D2EHPA as a carrier, due to its properties, was a 0.002 M solution of metal ions (initial pH = 4) without HCl used. As the stripping phase, a nitric acid solution with a concentration of 0.5 M and volume of 100 cm 3 was applied. The surface membrane area was equal to 17 cm 2 . The aqueous solutions were pumped from the external reservoirs by a peristaltic pump (GILSON MINIPULS 3) at a 16 cm 3 /min flow rate. During the transport experiments, 1 mL samples from the feed and the receiving solution were taken for analysis with the flame atomic absorption method using a SPECTRAA 20ABQ Varian spectrophotometer. The atomic absorption spectroscopy operating parameters are presented in the Supplementary File (Table S2). All the experiments were carried out at room temperature (25 ± 2 • C).

Membrane Preparation
PIMs were prepared according to the procedure described in [18,19] using a solution casting and solvent evaporation technique. Cellulose triacetate (Acros Organics, Morris Plains, NJ, USA), as a polymer matrix, and o-nitrophenyl octyl ether (Alfa Aesar, Kandel, Germany, 98%), as a plasticizer, were used. The composition and thickness of the membranes applied in the experiments are presented in Table 1.

Model Calculations
The ordinary differential equations (ODEs) were solved with the Berkeley Madonna program v.8.1 (Berkeley, CA, USA) using the Rosenbrock (stiff) method. This method is related to the Runge-Kutta method; however, it possesses excellent stability properties, is computationally efficient, and preserves the positivity of the solutions [38,39]. The parameter values were estimated using the curve fit procedure which minimizes the deviation between the model output and dataset. The dependences of concentration changes in the feed and stripping solution were fitted simultaneously.
The most common parameter used to select the model with the best fit quality is the determination coefficient. However, it was proven that for nonlinear models, the application of the determination coefficient for the model selection is questionable because R-squared does not distinguish between good and bad nonlinear models [40].
Therefore, as the fit quality parameter, the standard error of the regression (s y ) was calculated using the following equation: where RSS denotes the residual sum of squares, MSE is the mean squared error, y i is the observed value of the response variable,ŷ i is a predicted value of the response variable, n is the number of observations (the sample size), and p is the (total) number of estimated parameters. Moreover, the criteria for model selection among a finite set of models, such as the Akaike information criterion (AIC), Bayesian (Schwarz) information criterion (BIC), and Hannan-Quinn information criterion (HQC), were also calculated. These criteria are most frequently applied as measures of the goodness of fit of a statistical model and are defined by [40] where ln(L) denotes the log-likelihood function of the statistical model defined by The Hannan-Quinn information criterion was calculated using the following equation [41]: Note that the model with lower s y , AIC, BIC, and HQC values is preferred. Before the model selection, another test should be performed to examine the random distribution of the residuals. For this purpose, a nonparametric runs test was used [42]. The best model was therefore selected from those in which the randomness of the residuals was fulfilled.

Influence of the Feed and Stripping Solution Volume Changes
During the operation of the membrane system, samples of the feed and receiving solution were taken for concentration analysis (without return). Consequently, the effect of the change in the volume of the respective solutions was taken into account in the model calculations. The actual change in volume can be described by a step function, which is problematic to include in numerical calculations carried out in the Berkeley Madonna ODE solver. Therefore, a nonlinear function (third-degree polynomial) was used to describe the continuous change in volume during the system operation. As an example, in Figure 1a, the results of the calculations for the model proposed by Danesi are compared with those of the same model taking into account the influence of the feed and receiving solution volume (step and continuous). A comparison of the fits of the models represented by Equations (7) and (8) is presented in Figure 1b. distribution of the residuals. For this purpose, a nonparametric runs test was used [42]. The best model was therefore selected from those in which the randomness of the residuals was fulfilled.

Influence of the Feed and Stripping Solution Volume Changes
During the operation of the membrane system, samples of the feed and receiving solution were taken for concentration analysis (without return). Consequently, the effect of the change in the volume of the respective solutions was taken into account in the model calculations. The actual change in volume can be described by a step function, which is problematic to include in numerical calculations carried out in the Berkeley Madonna ODE solver. Therefore, a nonlinear function (third-degree polynomial) was used to describe the continuous change in volume during the system operation. As an example, in Figure 1a, the results of the calculations for the model proposed by Danesi are compared with those of the same model taking into account the influence of the feed and receiving solution volume (step and continuous). A comparison of the fits of the models represented by Equations (7) and (8) is presented in Figure 1b. The results of the concentration vs. time dependence calculations for the feed solution are almost the same regardless of whether the change in volume is taken into account or not. For the stripping solution dependence, the application of the continuous (solid line) or step function (dotted line) for the description of the volume changes leads to a much better fit of the model to the experimental data. The percentage error between the calculated values of concentration in the receiving solution for the two best-fit models does not exceed 1%. Therefore, the solution volume changes in the membrane system were described by the third-degree polynomial in further model calculations. An example of a third-degree polynomial fitting to the experimental results is presented in the Supplementary File ( Figure S1).

TOPO as a Carrier
The experimental results of Cd(II), Zn(II), Pb(II), and Cu(II) ion transport through PIMs with TOPO as a carrier were used for the best model selection. No transport of Cu(II) ions was observed in this study, i.e., the concentration was below the detection limit of AAS. The values of the fit quality parameters and the initial maximum fluxes are listed in Table 2. The lowest values are in bold. Table 2. The calculated fit quality parameters and initial maximum fluxes for the system with TOPO as a carrier.

Ion
No.
Scheme The results indicate that in the case of Cd(II) and Zn(II) ion transport, the lowest values of the s y , AIC, BIC, and HQC parameters were obtained for model No. 4, while model No. 3 was best-fitted to the experimental data of Pb(II) ion transport. Despite the similarity of the values of the initial fluxes calculated using different models, it should be remembered that the fit quality of the models is different. For example, Figure 2   The Zn(II) concentration differences in the feed solution, calculated by models No. 1 and No. 4, were small. Therefore, the same values of the initial maximum fluxes (JM = 3.400 × 10 −10 mol/cm 2 ·s) were evaluated. Significant differences are visible only in the concentration dependences for the receiving solution, leading to a maximum percentage error ex-  The Zn(II) concentration differences in the feed solution, calculated by models No. 1 and No. 4, were small. Therefore, the same values of the initial maximum fluxes (J M = 3.400 × 10 −10 mol/cm 2 ·s) were evaluated. Significant differences are visible only in the concentration dependences for the receiving solution, leading to a maximum percentage error exceeding 5% at the end of the transport. Such a value is inacceptable when the model is used in analytical applications, e.g., to predict the enrichment factor of analytes.

Aliquat 336 as a Carrier
The kinetic models were applied for the description of the transport of Cd(II), Zn(II), Pb(II), and Cu(II) ions through PIMs containing Aliquat 336 as a carrier. Similar to the system with TOPO as a carrier, transport of Cu(II) ions was not observed (below the detection limit of AAS). The calculated model fit quality parameters are presented in Table 3, along with the initial maximum fluxes. Based on these quantities, the most appropriate models were selected.  3 showed the best fit in the case of Cd(II) ion transport. In this membrane system, the best-fitted models led to initial maximum fluxes higher than 7 to 11%, especially for preferentially transported Cd(II) and Zn(II) ions.
A comparison of the fits of models No. 3 and No. 1 to the experimental data for Cd(II) ion transport is shown in Figure 3.
The experimental results and the fitted model No. 3 indicate the occurrence of a time lag in the investigated system and thus the accumulation of transported Cd(II) in the membrane. The same effect was also observed for Zn(II) and Pb(II) ion transport. This phenomenon is fundamental in analytical applications since accumulation reduces the transport efficiency and enrichment factor. Note that for models No. 1 and No. 2, the accumulation of the transported substance in the membrane was not taken into account.

Cyphos IL 101 as a Carrier
The fit quality measures obtained from fitting the kinetic models to the experimental data of Cd(II), Zn(II), Pb(II), and Cu(II) ion transport through PIMs containing Cyphos IL 101 as a carrier are presented in Table 4. The experimental results and the fitted model No. 3 indicate the occurrence of a time lag in the investigated system and thus the accumulation of transported Cd(II) in the membrane. The same effect was also observed for Zn(II) and Pb(II) ion transport. This phenomenon is fundamental in analytical applications since accumulation reduces the transport efficiency and enrichment factor. Note that for models No. 1 and No. 2, the accumulation of the transported substance in the membrane was not taken into account.

Cyphos IL 101 as a Carrier
The fit quality measures obtained from fitting the kinetic models to the experimental data of Cd(II), Zn(II), Pb(II), and Cu(II) ion transport through PIMs containing Cyphos IL 101 as a carrier are presented in Table 4.   The lowest values-which indicate a better fit of the model to the experimental data-were observed for model No. 3 used for the description of the Cd(II) and Pb(II) concentration vs. time dependences. In the case of Zn(II) and Cu(II) ion transport, only one model satisfied the criterion of the randomness of the data. For the Cu(II) ions, it was model No. 2, and for the Zn(II) ions, it was model No. 4. The results show that, similar to the system with Aliquat 336 as the carrier, for Cd(II), Zn(II), and Pb(II) ion transport, accumulation in the membrane was also observed. The biggest difference between the initial maximum flux values calculated from the appropriate models was found for the transport of Cu(II) ions; the J M calculated by the correctly fitted model was over 55% higher than the values calculated by the other models. The fitting of the worst (model No. 1) and best models (No. 4) to the experimental data of Zn(II) ion transport is compared in Figure 4. one model satisfied the criterion of the randomness of the data. For the Cu(II) ions, it was model No. 2, and for the Zn(II) ions, it was model No. 4. The results show that, similar to the system with Aliquat 336 as the carrier, for Cd(II), Zn(II), and Pb(II) ion transport, accumulation in the membrane was also observed. The biggest difference between the initial maximum flux values calculated from the appropriate models was found for the transport of Cu(II) ions; the JM calculated by the correctly fitted model was over 55% higher than the values calculated by the other models. The fitting of the worst (model No. 1) and best models (No. 4) to the experimental data of Zn(II) ion transport is compared in Figure 4.

D2EHPA as a Carrier
The calculated fit quality characteristics of the analyzed models used for the description of Cd(II), Zn(II), Pb(II), and Cu(II) ion transport through PIMs containing D2EHPA as a carrier are compared in Table 5. Table 5. The calculated fit quality parameters and initial maximum fluxes for the system with D2EHPA as a carrier.

Ion
No.
Scheme  The transport kinetics in this system can be described by various models satisfying the criterion of the randomness of the data, especially for the transport of Cd(II) ions. From the comparison of the fit quality parameter results, all of the fit criteria using model No. 3 were lower than those of the other models for Cd(II) and Zn(II) ion transport. The lowest s y , AIC, BIC, and HQC values for Pb(II) ion transport were found for model No. 4, while for Cu(II) ions, the lowest values were found for model No. 2. The calculated initial maximum fluxes using different models were similar. However, a decidedly different fit of the models to the experimental data was observed (see Figure 5). For Cu(II) ion transport, the initial maximum flux calculated by model No. 2 was approx. 14% higher when compared to the values calculated by the other models.

Reactive Ionic Liquid (RILC8_Br) as a Carrier
A summary of the fit quality measures of the analyzed models is presented in Table  6. In this membrane system, no transport of Cu(II) ions was observed.

Reactive Ionic Liquid (RILC8_Br) as a Carrier
A summary of the fit quality measures of the analyzed models is presented in Table 6. In this membrane system, no transport of Cu(II) ions was observed.
All the fit quality measures indicate that there was only one kinetic model appropriate for the description of Cd(II), Zn(II), and Pb(II) ion transport through PIMs using RILC8_Br as a carrier. In the case of Cd(II) and Zn(II) ion transport, except for model No. 2, no other model met the random distribution of residuals criterion. Choosing an inappropriate model leads to substantial errors in the initial maximum flux values. The calculated J M values from the best-fit model were higher by approx. 22% for Cd(II) ions and up to 74% for Zn(II) ions. The proper choice of this model for the description of metal ion transport through PIMs containing reactive ionic liquids was also confirmed by the results presented in [18]. The fit of models No. 2 and No. 4 to the experimental data of Cd(II) ion transport is shown in Figure 6. The results indicate a much better fit of model No. 2. Moreover, no time lag effect was observed in the investigated system. The undetectable time lag may result primarily from the properties of the carrier used, as well as the application of the thinnest membrane in the experiments. transport through PIMs containing reactive ionic liquids was also confirmed by the results presented in [18]. The fit of models No. 2 and No. 4 to the experimental data of Cd(II) ion transport is shown in Figure 6. The results indicate a much better fit of model No. 2. Moreover, no time lag effect was observed in the investigated system. The undetectable time lag may result primarily from the properties of the carrier used, as well as the application of the thinnest membrane in the experiments.

Conclusions
An appropriate mathematical model enables a quantitative prediction of concentration changes in external solutions and the calculation of, e.g., the time needed to achieve the maximum concentration, recovery, or enrichment factor of the transported substances.
Among the numerous equations used to describe transport kinetics, it is impossible to choose only one as the most appropriate. The presented results indicated that the application of the most frequently used model (proposed by Danesi) is significantly limited because of the nonrandomly distributed residuals.
The selection of the appropriate model should be based on the values of the standard error of the regression, Akaike information criterion (AIC), Bayesian (Schwarz) infor-

Conclusions
An appropriate mathematical model enables a quantitative prediction of concentration changes in external solutions and the calculation of, e.g., the time needed to achieve the maximum concentration, recovery, or enrichment factor of the transported substances.
Among the numerous equations used to describe transport kinetics, it is impossible to choose only one as the most appropriate. The presented results indicated that the application of the most frequently used model (proposed by Danesi) is significantly limited because of the nonrandomly distributed residuals.
The selection of the appropriate model should be based on the values of the standard error of the regression, Akaike information criterion (AIC), Bayesian (Schwarz) information criterion (BIC), or Hannan-Quinn information criterion (HQC) after the runs test evaluation (residuals' randomness check).
Of the models presented in this report, model No. 4 is the most universal. However, the model selection should be individualized for each experimental relationship. It was also found that a nonlinear equation (third-degree polynomial) can be successfully used to describe solutions' volume changes in a membrane system (because of sampling), leading to a better fit of the model to the experimental data.
The results indicated that the models that have not been used thus far for transport description in PIMs, i.e., models No. 3 and No. 4, can be successfully applied. These models are particularly important in the case of systems where a time lag is observed. This mainly applies to systems with carriers characterized by a high partition coefficient (high sorption of transported substances to the membrane), slow diffusion inside the membrane (e.g., because of the relatively high membrane thickness), or slow kinetics of extraction and re-extraction at the respective membrane interfaces.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/membranes13020236/s1, Table S1. List of differential equations used in an appropriate models. Table S2. Atomic absorption spectroscopy operating parameters. Table S3. The calculated kinetic parameters and initial maximum fluxes for the system with TOPO as a carrier. Table S4. The calculated kinetic parameters and initial maximum fluxes for the system with Aliquat 336 as a carrier. Table S5. The calculated kinetic parameters and initial maximum fluxes for the system with Cyphos IL 101 as a carrier. Table S6. The calculated kinetic parameters and initial maximum fluxes for the system with D2EHPA as a carrier. Table S7. The calculated kinetic parameters and initial maximum fluxes for the system with reactive ionic liquid (RILC8_Br) as a carrier. Figure S1. The 3rd degree polynomial fitting to the experimental results for system with D2EHPA as a carrier. Figure S2. The best fitted models for predicting Cd(II) (a), Zn(II) (b), and Pb(II) (c) ions transport through PIM with TOPO as a carrier. Figure S3. The best fitted models for predicting Cd(II) (a), Zn(II) (b), and Pb(II) (c) ions transport through PIM with Aliquat 336 as a carrier. Figure S4. The best fitted models for predicting Cd(II) (a), Zn(II) (b), Pb(II) (c), and Cu(II) (d) ions transport through PIM with Cyphos IL 101 as a carrier. Figure S5. The best fitted models for predicting Cd(II) (a), Zn(II) (b), Pb(II) (c), and Cu(II) (d) ions transport through PIM with D2EHPA as a carrier. Figure S6. The best fitted models for predicting Cd(II) (a), Zn(II) (b), and Pb(II) (c) ions transport through PIM with a reactive ionic liquid (RILC8_Br) as a carrier.
Funding: This research received no external funding.